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Abstract 

In this paper we investigate the computational aspects of some recently proposed 
iterative methods for approximating the canonical tight and canonical dual window 
of a Gabor frame (g,a,b). The iterations start with the window g while the itera- 
tion steps comprise the window g, the k th iterand 7^, the frame operators S and 
Sk corresponding to (g,a,b) and (7^,(2, 6), respectively, and a number of scalars. 
The structure of the iteration step of the method is determined by the envisaged 
convergence order m of the method. We consider two strategies for scaling the terms 
in the iteration step: norm scaling, where in each step the windows are normalized, 
and initial scaling where we only scale in the very beginning. Norm scaling leads to 
fast, but conditionally convergent methods, while initial scaling leads to uncondi- 
tionally convergent methods, but with possibly suboptimal convergence constants. 
The iterations, initially formulated for time-continuous Gabor systems, are con- 
sidered and tested in a discrete setting in which one passes to the appropriately 
sampled-and-periodized windows and frame operators. Furthermore, they are com- 
pared with respect to accuracy and efficiency with other methods to approximate 
canonical windows associated with Gabor frames. 
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1 Introduction 

We consider in this paper iterative schemes for the approximation of the canon- 
ical tight and canonical dual windows associated with a Gabor frame. We 
refer to [6, Ch. 8-10] and [11, Ch. 5-9, 11-13] for recent and comprehensive 
treatments of the theory of Gabor systems and frames; to fix notations and 
conventions we briefly give here the main features. We denote for g G L 2 (R) 
and a > 0, b > by (g, a, b) the collection of time-frequency shifted windows 

g n a,mb, m,neZ, (1.1) 

where for x, y G R we denote 

9x, y = e 2aiyt g(t-x), teR. (1.2) 

We refer to g as the window and to a and b as the time-shift and frequency- 
shift parameters, respectively, of the Gabor system (g,a,b). When there are 
A > 0, B < oo, called the lower and upper frame bound, respectively, such 
that for all f £ L 2 (R) there holds 

oo 

^ll/H 2 < E K/,^m6)| 2 <S||/|| 2 , (1-3) 

m,n=— oo 

we call (g,a,b) a Gabor frame. When in (1.3) the second inequality holds for 
all / G L 2 (R), we have that 

oo 

/ G L 2 (R) ^ Sf := ]T (f,g na , mb )g na , mb (1.4) 

m,n=— oo 

is well-defined as a bounded, positive, semi-definite linear operator of L 2 (R). 
We call S the frame operator of ((?, a, 6). The frame operator commutes with 
all relevant shift operators, i.e., we have for all / G L 2 (R) 

Sfna,rnb = (Sf) namb , m, G Z. (1.5) 

We shall assume in the remainder of this paper that (g, a, b) is a Gabor frame. 
Thus the frame operator S is positive definite and therefore boundedly invert- 
ible. There are two windows canonically associated to the Gabor frame (g, a, b). 
These are the canonical tight window g t and the canonical dual window g d , 
defined by 

g t = S- 1 ' 2 g, S i d = S- 1 g, (1.6) 
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respectively. The practical relevance of these windows is that they give rise to 
Gabor series representations of arbitrary /6i 2 (R) according to 



oo oo 
/ y ' (/ ' 9na,mb) 9na,mb ^ ' (/ ' 9na,mbj 9na,mbi (1.7) 



m,n=— oo 



where both series are L 2 (M)-convergent. Furthermore, the Gabor systems 
(<7*, a, b) and (^g d , a, bj are Gabor frames themselves with frame operators equal 
to the identity / and S' 1 , respectively. 

The computation of g l and g d according to (1.6) requires taking the inverse 
square root and the inverse of the frame operator S, respectively. In the often 
occurring practical case that ab is a rational number, the frame operator is 
highly structured which allows relatively efficient methods for computing S* -1 , 
see [27]. The computation of S~* is much more awkward, even in the case that 
ab is rational, and often requires advanced techniques from numerical linear 
algebra, see for instance [13]. 

In [20] the calculus of Gabor frame operators was combined with a use of 
the spectral mapping theorem and Kantorovich's inequality to analyze an 
iteration scheme for the approximation of g l that was proposed around 1995 
by Feichtinger and Strohmer (independently of one another). In this iteration 
scheme one sets 7 = g and for k — 0, 1, . . . 



L lk+1 = \ctklk + \ faS k Sfc, «fc = tt-t v Pk = — — — • ( i ■< 



— re IK ' r-rc~fc IKl —re II || 7 /-re i 

2 2 || 7fc || s k ^ k 

where S k is the frame operator corresponding to (7^,0,6). It was shown in 
[20] that (7^, a, b) is indeed a Gabor frame, and that tP^m converges (at least) 

quadratically to From the numerical results in [20] for the case that g is 

the standard Gaussian window 2 1//4 exp (— irt 2 ) and a = b = l/y/2 it appears 
that the resulting method compares favorably with other iterative techniques 
for computing inverse square roots, [24,25,29]. 

The investigations in [20] were followed by the introduction in [18,19] of two 
families of iterative algorithms for the approximation of g l and g d , in which the 
iteration step involves the initial window g and the frame operator S as well 
as the current window 7^ and frame operator S k , but frame operator inversion 
as in (1.8) do not occur. The following instances of these two families were 
analyzed in [18,19]. Again we set 70 = g, and for k — 0, 1, . . . 



3 1 11 

! ' lk+i = -®klk - -PkSklk] OLk = -fi — m, Pk = Yd iT' ( L9 ) 

2 2 llTfcH ||<Vyfc|| 
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_ 15 5 3 c2 

lk+1 — —ZkOlk — -^kl^klk + ~^k2^klk] 
11 1 

£ko — 71 — [7, £fci — Yo FT' £fc2 ~~ TTc2 — m' (1-10) 

llTfcll IWfcll ||5fc7fc|| 

for the approximation of and 



V. 



7fc+i = 34o7fc - 35 k iS k g + 5 k2 SS k ^ k ; 

4o = 7i — [7, 4i = 7775 — n") ^2 = 777575 — if) (1-12) 
llTfcll l™ll \\SSklk\\ 



for the approximation of g d . 

The algorithms II- V are, contrary to algorithm I, conditionally convergent in 
the sense that the frame bound ratio of the initial Gabor frame (g, a, b) 
should exceed a certain lower bound. Accordingly, in algorithm II and III we 
have that converges to quadratically and cubically when > \ and 

^ > f , respectively. In algorithm IV and V we have that ypj^j converges to jj^yj 

quadratically and cubically when > \ (V§ — l) and I > 0.513829766 . . ., 
respectively. A remarkable phenomenon that emerged from the preliminary 
experiments done with the algorithms around 2002, was the fact that the 
lower bounds for the algorithms II, III seem far too pessimistic while those for 
the algorithms IV and V appear to be realistic. 

In the algorithms just presented, all computed windows are normalized. We 
shall refer to this as norm scaling. Another possibility that we will investigate, 
is to replace all scalars ce's, /3's, e's and 5's that occur in (1.9-1.12) by 1, and 
then initially scale the windows by replacing 



g by g/B 1 ' 2 ,S by S/B. (1.13) 

We shall refer to this scaling strategy as initial scaling. If B is (an estimate 
for) the best upper frame bound maxcr(S') then the algorithms will be un- 
conditionally convergent, with guaranteed desired convergence order, but with 
convergence constants that may not be as good as the ones that can be ob- 
tained by using the norm scaling as described by (1.9-1.12). 

Matrix versions of algorithms 1 1- 1 1 1 without any scaling have been treated in 
[5]. In [14,23] the matrix version of algorithm II is considered using norm scal- 
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ing and a scaling method that approximates the optimal scaling. The matrix 
version of algorithm IV is known as a Schulz iteration, see [28]. The fact that 
S, and therefore ip (S) with ip continuous and positive on the spectrum of S, 
commutes with all relevant shift operators, allows us to formulate the iteration 
steps on the level of the windows themselves. 

In this paper we investigate the algorithms II- V, using both norm and initial 
scaling, with more emphasis on computational aspects than in [20,18,19]. Here 
it is necessary to consider sampled-and-periodized Gabor systems in the style 
of [16]. This allows for a formulation and analysis of the algorithms I-V in 
an entirely similar way as was done in [20,18,19]. Thanks to the fact that the 
involved (canonical) windows and frame operators behave so conveniently un- 
der the operations of sampling and periodization, the observations done on the 
sampled-and-periodized systems are directly relevant to the time-continuous 
systems. We must restrict here to rational values of ab, and this gives the frame 
operator additional structure which can be exploited in the computations as 
dictated by the recursion steps, also see [2,27,30] for this matter. 

The notions "smart" (or, rather, "smart but risky") and "safe" (or, rather, "safe 
but conservative") were introduced in a casual way in [19] to distinguish be- 
tween cases where the stationary point (s) of the function transforming (frame) 
operators according to (4.2), (4.9) has a good chance to be well-placed in the 
middle of and on the "safe" side of the relevant spectral set, respectively. In 
the present paper, we choose to refer to the strategies leading to smart and 
safe modes as "norm scaling" and "initial scaling", respectively, and discard 
the terms "smart" and "safe" altogether. 



2 Paper outline and results 

In Section 3 and 4 we present the basic results of [18,19] on transforming g 
into 7 = ip (S) g, where <p is a function positive and continuous on a (S), so as 
to obtain a window 7 whose frame operator S 7 (for approximating g v ) or the 
operator (SSj) 1 ^ 2 =: Z 7 (for approximating g d ) is closer to (a multiple of) the 
identity / than S itself. Here we recall that (</*, a, b) has frame operator / and 
that (^g d ,a,bj has frame operator S* -1 . Thus we relate the operators S and 
S 7 and their frame bounds, and we present a bound for the distance between 
(the normed) 7 and g l in terms of the frame bounds of S 1 . Similarly, we relate 
the frame bounds of (g, a, b) and the minimum and maximum of a (Z y ), and 
we present a bound for the difference between (the normed) 7 and g d in terms 
of the latter minimum and maximum. Next, in Section 4, the choice of <p is 
specified so as to accommodate the recursions of type II, III and of type IV, 
V which gives us a means to monitor the frame bound ratio (for g l ) and 
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of the ratio between minimum and maximum of a (Z^) (for g d ) during the 
iteration process. 



In Section 5 we present the algorithms using only initial scaling. 



In Section 6 we elaborate on the observation that all algorithms take place 
in the closed linear span C g of the adjoint orbit [gj/b,i/ a j,l £ ^j. Here the 
dual lattice representation of frame operators is relevant as well as an operator 
norm to measure the distance of Sk (for g l ) and of (SSk) 1 ^ 2 (for g d ) from (a 
multiple of) the identity. The consideration of the algorithms in the space L g 
reveals a fundamental difference between the algorithms for computing g l and 
g d that manifests itself in the totally different after-convergence behaviour of 
the two families of algorithms. 



In Section 7 we give some considerations in the Zak transform domain, so as 
to produce examples of Gabor frames for which a specific algorithm diverges. 



In Section 8 we discuss the discretization and finitization aspects (through 
sampling and periodization) that have to be taken into account since the 
algorithms are to be tested numerically. 



In Section 9 we show how the algorithms can be expressed for discrete, finite 
Gabor systems, and show that the algorithms are scalar iterations of the sin- 
gular values of certain matrices. We present an efficient implementation of the 
iterative algorithms, and we list the window functions we have used to test 
the algorithms. 



In Section 10 we present our experimental results, compare them with what the 
theory predicts and with other methods to compute tight and dual windows. 
We provide examples that show the quadratic and cubic convergence of the 
algorithms, and the exponential divergence of the dual iterations after the 
initial convergence. We give an example that breaks the norm scaling schemes 
for both the tight and dual iterations, and show how various error norms of 
the iteration step behave. Comparisons with other methods are made: We 
show that the tight iterations are competitive with respect to computing time 
and superior with respect to precision. Finally, we show that the number 
of iterations needed for full convergence of the algorithms are dependent on 
the frame bound ratio, but independent of the structural properties of the 
discretization. For initial scaling, we show that it is easy to choose a scaling 
parameter that gives almost optimal convergence. 
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3 Frame operator calculus and basic inequalities 



The basic theory to analyze the recursions appears somewhat scattered in 
[20,18,19]; for the reader's convenience, we give in Section 3 and 4 a concise 
yet comprehensive summary of the basic results and ideas. We let (g, a, b) be a 
Gabor frame with frame operator S and best frame bounds A = min a (S) > 0, 
B = max a (S), where a (S) is the spectrum of S. In this section we present the 
basic inequalities expressing the approximation errors in terms of the (frame) 
bounds on the involved (frame) operators. These inequalities are a consequence 
of the calculus of Gabor frame operators, the spectral mapping theorem and 
Kantorovich's inequality. 

Proposition 1 Let ip be continuous and positive on [A,B], and set 7 := 
ip (S) g. The following holds. 

(i) (7, a, b) is a Gabor frame with frame operator S 1 := Sip 2 (S) and best 
frame bounds 



A„ := min sip 2 (s) , B~ := max sip 2 (s) . 

1 s£a(S) 1 s£a(S) 



(3.1) 



Furthermore, 
and 
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f = S-^g = SZ^ = 



9 



1711 

(ii) Let Z 7 := (SS 7 ) 1/2 = Sip {S), and 



O =^ 

-D-y 



(3.2) 
(3.3) 



E 1 := min a (Z 7 ) = min sip(s), 

sEa(S) 

F 7 := max a (Z 7 ) = max sip (s) . 

setr(S) 



(3.4) 
(3.5) 



Then 



S , ^ 7 7 



S-yg, 



(3.6) 



and 



1_ 
\l\ 



g u 



< (1 - < 2 ) 



£ty — — . 



(3.7) 



The basic result (i) gives us a clue how to produce a good approximation tAt 

of jj^jj-: Take ip such that sip 2 (s) is flat on a (S) C [A, B] so that the number 
Q 1 in (3.3) is close to 1. Similarly, by the basic result (ii), the number i? 7 in 
(3.7) is close to 1 when ip is such that sip (s) is flat on a (S) C [A, B], and 
then we obtain a good approximation of -tt^tt. In the next two subsections, we 
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use this basic result repeatedly with polynomials cp of fixed degree m so as to 
obtain iterative approximations of g f and g d . 



4 Norm scaling 



4-1 Iterations for approximating g t 



We consider iteration schemes 



lo = g; lk+i = ¥k{S k )lk, A; = 0,1,..., (4.1) 

for the approximation of (?*, where Sk is the frame operator of (7^,0,6). We 
use here the basic result (i) repeatedly with 



# = 7fc, S = S k and 7 = j k+1 , S 7 = S k+1 = S k <p 2 k (S k ) . (4.2) 
For k — 0, 1, . . . we have that 7^ = and that 

where and are the best frame bounds of Sk- The numbers A k , Bk can 
be computed and estimated recursively according to A = A, B = B and 



Ik _ J_ 



A k +i = mm sip 2 k (s) > min sip 2 k (s) , 

sea(S k ) s£[A k ,B k ] 

B k+ i = max s^fc (s) < max sy?); (s) , 

sea(S k ) s£[A k ,B k ] 

for fc = 0, 1, . . .. 



(4.4) 
(4.5) 



We should choose ip k such that sip k (s) is flat on a (S k ). To that end there is 
proposed in [19, Subsec. 5.1] for m — 2, 3, . . . the choice 



m—1 



¥k (s) = amjakjS 3 , a kj = S 3 k ^ k 
3=0 



where the defined by 



(4.6) 




\{l-x) =2^ a mj x J , x>0. 
I / 3=0 



(4.7) 
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The motivation for this choice is as follows. The left-hand side of (4.7) is the 
(m — l) th order Taylor approximation of x^ 1 ^ 2 around x — 1, while (when Q k 
is sufficiently close to 1) 

Hence s<p k (s) = (s 1 / 2 ^ (s)) should be expected to be flat on a(S k ), with 



1 — Qk+i potentially of order (1 — Q k ) m ■ 

When m = 2, 3 we get the iterations II, III in (1.9), (1.10). It is shown in [18, 
Sec. 4] and [19, Sec. 6] that for m = 2 the quantity Q k = ^ increases to 1 and 

that jj^jy — * quadratically when k — > oo, provided that ^ > |. For m = 3 

it is shown in [18, Sec. 8] that Q k increases to 1 and that — > cubically 

when /c — > oo, provided that > |. In [18,19] it was observed for m = 2, 3 
that the choice of a k j causes sip 2 , (s) to have one or more stationary points in 
[Ak, Bk] so that the odds for flatness of s(p k (s) on [A k , B k \ are favourable. 

4-2 Iterations for approximating g d 

We consider iteration schemes 

7o = g\ ik+i = y?fc (Z k ) ik, k = 0, 1, . . . , (4.9) 

for the approximation of g d , where Z k = (SSk) 1 ^ 2 with S k the frame operator 
of (7fe,a, b). It is seen from (4.9) by induction that ^ k = vp k (S)g for some 
function ip k . Hence by the basic result (i) with <p = tp k , we have that 

S k = Sil> 2 k (S),Z k = Srl> k (S), (4.10) 
and, by the basic result (ii), that 



Ik 9 d 



hk\\ \\g d \ 



^-*n^ R *=% (4 - n) 



where E k = miner (Z k ), F k = maxa (Z k ). A further use of the calculus of 
frame operators as given by the basic result (i), yields 

Z k+ i = Z k ip k (Z k ) . (4.12) 

Consequently, the numbers E k , F k can be computed and estimated recursively 
according to E = A, F = B and 



9 



E k+1 = min zip k (z) > min zip k (z) , 

zea(Z k ) ze[E k ,F k ] 

F k+ i = max z<p k (z) < max zip k (z) , 

zea(Z k ) ze[E k ,F k ] 



(4.13) 
(4.14) 



for k = 0,1,.... 

We should choose (p k such that z(p k (z) is flat on a (Z k ). To that end there is 
proposed in [19, Subsec. 5.2] for m — 2, 3, . . . the choice 



m— 1 

V?fc (*) = H b mj (5 kj z\ 

3=0 



Pkj — Z k lk 



(4.15) 



where the 6 m ,- are defined by 



m— 1 



m— 1 



(! - X ) Z = S bmjX 3 , x>0. 

1=0 3=0 



(4.16) 



The motivation for the proposal is similar to the one for the choice of (p k in (4.6) 
in Subsec. 4.1; we now note that the left-hand side of (4.16) is the (m — 1)* 
order Taylor approximation of x^ 1 around x — 1. The implementation of the 
resulting recurrence step 



m—i Z k ^ k i 2 

7fc+l = ^2 bmj J^~: Z k = (SSk) 1 ^ 2 , 

Z klk 



3=0 



(4.17) 



is made feasible by the observation that, thanks to the second item in (3.6), 
Z k -i k = S k g so that 



Zfc r 7fc = {SS k ) r r/ k , Zl T+1 ^ k = (SS k ) r S k g, r = 0,l, 



(4.18) 



When m = 2,3 we get the recurrences IV, V in (1.11) and (1.12). It is shown 
in [18, Sec. 5] and [19, Sec. 7], that for m = 2 the quantity R k = |£ increases 



to 1 and that 

\hk\ 



quadratically when k — > oo, provided that 



> 



| (y/E — l). For m = 3 it is shown in [18, Sec. 9] that R k increases to 1 and 



2 

that 



Ik 



cubically when k — > oo, provided that -= > 0.513829766 . 



In [18,19] it was observed for m = 2, 3 that the choice of (5 k j causes zip k (z) to 
have one or more stationary points in [E k , F k }. 
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5 Initial scaling 



The algorithms II- V are guaranteed to converge when the lower bound ratio ^ 
of (g, a, b) exceeds a certain value. The proofs, as given in [18] and [19], require 
a careful analysis of the extreme values of the functions <fk on the spectra of 
the relevant operators and can become quite complicated, especially in the 
cases of algorithms III, V. However, the algorithms are efficient in the sense 
that the envisaged convergence order to is realized with favourable convergence 
constants. In practice, as the experiments in Section 10 show, the algorithms 
II, III turn out to converge in almost all cases, even when the frame bound 
ratio is close to 0. However, divergence of the algorithms IV, V occurs much 
more frequently. In Section 7 we present examples, using the Zak transform, 
of frames (g, a, b) such that algorithm II and IV diverges. 

It would be desirable to have versions of the algorithms that are guaranteed 
to converge, no matter how small the frame bound ratio of the initial frame 
is (as long as it is positive). In the following, we present the initial scaling 
versions of the algorithms that converge at the envisaged convergence order 
to, possibly with suboptimal convergence constants. Since we can freely switch 
scaling strategy, a possible strategy is to initially scale such that convergence is 
guaranteed, and to continue until one is confident that the relevant condition 
number exceeds the specific lower bound so that the norm scaling mode can 
be applied from that point onwards. 

The introduction in [19] of the notion of "safe modes" was prompted by an ob- 
servation by M. Hampejs who prescaled the window g (and the frame operator) 
and deleted all normalization operations in the recursion step of algorithms II, 
IV. In the present paper, the prescaling is done in such a way that the scaled 
S has its spectrum exclusively in the attraction region of the function ip de- 
scribing the simplified recursion. More specifically, we consider the iteration 
steps as given in Subsections 4.1, 4.2, with all a's and /3's equal to 1. The </?'s 
thus obtained are independent of k and are given by 



m—1 
1=0 




m—1 



1 - S) = £ V-rajS 1 , S > 0, 
3=0 



(5.1) 



and 



m—1 m—1 



E (!" *) = EW. ^ >0 > ( 5 - 2 ) 

1=0 j=0 

respectively. The relevant spectra transform by the spectral mapping theorem 
according to 
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1 
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m=2 . 
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- 
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y i 


I s il I 



(a) s > ^ s (s)) 



2 2.5 

2 



0.5 1 1.5 2 2.5 3 



(b) z > ^ (z) 



Figure 1. The figure shows the two set of functions governing the convergence of the 
tight and dual iterations using initial scaling for order m = 1,2,3,4 with (p^ and 
tfm defined by (5.1) and (5.2). For (a), the attraction point 1 has attraction regions 
m = 2 : (0,3), m = 3 : (0,7/3), m = 4 : (0,2.525847988). For (b), the attraction 
point 1 has attraction region (0, 2) for m = 2, 3, 4. 



a(S k 



s e a(S k ) \ = a(S k+1 ) , 



(5.3) 



and 



zea(Z k )j=a(Z k+1 ), (5.4) 
respectively. 

The functions (p^ and tp d m are (m — l) th order Taylor approximations of s~ 1//2 
and z~ l around s = 1 and -2 = 1, respectively. Hence s (y?^ (s)) 2 and z^ d m (z) 
approximate 1 around s = 1 and z — 1, respectively. In Fig. 1 we have shown 
plots of the mappings 



a (Z k ) -> ^ (z) 



s>0^ s (^( s )) ,z>0^^(z) (5.5) 

for m = 1,2,3,4, respectively. Fig. 1(a) also appears in [5]. In all cases, the 
point s = 1 = z is an attractor for the region (0,2). Consequently, when 
S = S, Z = S have spectrum in (0,2), the spectra a(S k ), a (Z k ) converge 
to 1 as k — > oo, and the convergence is of order m in the sense that the 
ratio of minimum and maximum of the spectra converge to 1 at order m. 

1 /2 

Thus we should replace g by gj {B^ and S by S/B where B is such that 
a (S/B^J C (0, 2) to obtain iterations having m th order convergence to g l and 

/m1/2 , 

to [B) g , respectively. 
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Table 1 



Method. 


B or F 


I. 


VAB 




II. 


l(A + VAB + Bj 




III. 


±(B + A) + y±(Bi + A*) + ±(B 


-Af 


IV. 


\{E + F) 




V. 


\{F + E) + \yl\{F* + EP) + \{F- 


- E) 2 



The table shows the optimal scaling constant, B or F, for doing initial scaling of 
the five iteration types. 



To guarantee convergence an estimate of max a (S) is needed. In [2] a number 
of upper bounds of maxcr (S) are developed for discrete-time, periodic Gabor 
systems. A convenient upper bound for our purposes follows from the dual 
lattice representation of the frame operator S, see Section 6 for more details, 

as 

(5.6) 



max (j(S)<-rJ2\ {d, 9j/b,i/a) 
3,1 



We make some comments for scaling optimally in the first iteration step. We 
shall refer to this method as initial optimal scaling. Assume that a (S) consists 
of the entire interval [A, B\. Consider the tight iterations as described in this 
Section, and assume that we replace S by S/B. Then 



A 1 = min I s(y m (sj) 
B l = max I s (</4 ( 



s e 



s e 




(5.7) 
(5.8) 



Initial optimal scaling occurs for that value of B for which the ratio A\jB\ is 
maximal. The optimal value of the scaling parameter F for the dual iterations 
is defined in a similar way. For the five iteration types, the optimal B or F 
is shown in Table 1. All these numbers are close to the center of the interval 
[A, B] or [E,F\. It is not hard to show that all optimally scaled operators S 
and Z have their spectra in the attraction regions given in Figure 1 for the 
algorithms II- V. 

If we scale optimally in each iteration step, and not just the first, we get the 
best possible convergence constants. However, this method is not practically 
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feasible because of the repeated calculations of frame bounds, and we shall 
use it only as a reference method. We refer to it as constant optimal scaling, 
see Fig. 5. 



6 Considerations in the adjoint orbit space 



We consider the closed linear span C g of the adjoint orbit {gj/b,i/ a ]\l £ 
According to the duality principle of Gabor analysis we have that the adjoint 
orbit is a Riesz basis for C g (for this matter we refer to [6, Sees. 3.6 and 9.2], 
and [11, Ch. 7]). Furthermore, when / el 2 (M), the orthogonal projection of 
/ onto C g is given by 

v= s e (/• (A*J = s £ (** W) (A™. <"> 

As a consequence of 7^ = in all our algorithms, we see that 7^ G £ ff . There 
is also the Wexler-Raz biorthogonality relation, 

where 5 is Kronecker's delta. Finally, there is the following fundamental iden- 
tity of Gabor analysis. Assume that /, £, 7, h e (R) and that the three Ga- 
bor systems (/, a, 6), (£, a, 6), (7,0,6) have finite upper frame bounds. Then 
we have 

7j/6,J/a) (/j/6,I/o, (6-3) 

with absolute convergence at either side. We can regard (6.3) as a representa- 
tion result for the frame-type operator 



S-y,Z '■ f ~* ^2 (/' 7na,mf>) £na,mf>, (6.4) 



viz. as 



5 7,£ = i S 7j/6,i/a) (6-5) 



where [/^ is the unitary operator 



14 



U 3,l : f fj/b,l/a- (6.6) 

This is the dual lattice representation (also known as the Janssen representa- 
tion, see [6, Sec. 7.2] and [11, Corr. 9.3.7]) of the frame operator. In order for 
(6.5) to be well-defined, we assume that 7, £ satisfies the so-called Condition 
A': 

A': 

E|(f,7j/6,!/o)| < 00, (6.7) 
3,1 

see [11, Def. 7.2.1]. If £ = 7 then this is the Condition A introduced by 
Tolimieri and Orr in [32] . We refer to Appendix A where an instance, relevant 
in the present context, of a pair £,7 satisfying condition A' is given. 



6. 1 Estimate for upper frame bound 



If g satisfies Condition A, the frame operator S of (g, a, b) has the representa- 
tion 



S = ^ Yl {9, 9j/b,i/a) U jt i, (6.8) 



3,1 



with absolute convergence in the operator norm. Therefore, there is the upper 
bound 



^ = 71 E I {9, 9j/b,i/a) I (6-9) 
3,1 

for the best upper frame bound max a (S) of (g, a, b). 



6.2 Error measure 



We measure convergence of 7& to g l and g d by inspecting L 2 -distances of the 
normed windows. This quantity is bounded in terms of the numbers Qk and 
Rk in (3.3) and (3.7) that measure how close the operators Sk and Zk are to 
being a multiple of the identity operator. In the converse direction, it would 
be useful to have a measure on the windows that translates directly to the 
distance of Sk and Z k to (a multiple of) the identity operator. Such a measure 
can indeed be found. As to g t we note that when g satisfies Condition A then 
Sk has the representation 
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s * = ^ 12 (cik, (ik) j/bi i /a ) u jt i, 

3 J 



(6.10) 



whence 



5*-^ INI 2 / 



|(7fe>(7fc),/ M/a ) 



(6.11) 



As to g d we note that Z fe = S^fe (S*) and with ^ k = ip k (S) g there holds by 
frame operator calculus 



S<p k (S) / = £ (/» (7*) na,mb J 9na,mb- 

m,n 

Hence there is the representation 



(6.12) 



Z k = S<p k (S) = —J2 (q, (lk) j/b , l/a ) Uj,i- 
3,1 



Therefore 



Zk ~ ab <K9 ' lk ^ >I 



<l~ h E |(ff>(7*) j/6lI /a) 



(6.13) 



(6.14) 



Note that the quantities of the right-hand sides of (6.11) and (6.14) measure 
to what extent the Wexler-Raz condition (6.2) is violated. In Appendix A it is 
shown for a G N, b^ 1 G N and g satisfying condition A that the 7^ occurring in 
(6.11) and the g, ^ k occurring in (6.14) satisfy condition A and A', respectively. 
We shall refer to the right hand sides of (6.11) and (6.14) as the dual lattice 
norm. 



6.3 Influence of out- of- space components 



We have seen that all iterands 7^ of the algorithms are in C g . We briefly 
comment on the impact on the algorithms of 7^ having non-zero components 
orthogonal to C g (one can think here of round-off errors generating these com- 
ponents). To that end we consider the algorithms II and IV (assuming appro- 
priate scaling has been carried out), and we assume that they have converged 
to the extent that the operators S k and Z k agree within machine precision 
with the identity operator. 



As for algorithm II we thus have that 



3 l q 

7fc+i = o7fe - o^fc7fc = Ik 



(6.15) 
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within machine precision. Hence, possible out-of-space components in ^ k are 
reproduced within machine precision. As a consequence, we should expect that 
the error stays at its converged level when the iteration is continued beyond 
the point where machine precision is reached. 

Next we consider algorithm IV using initial scaling so that 

7 fc+ i = 27 fc - S k g. (6.16) 
The term Skg has the representation 



S *9 = ^ E (7fc, (lk) j/b , l/a ) 9j/b,i/a e C g . (6.17) 

3,1 

Furthermore, 



P^^E^^^Jg^/a. (6.18) 

3,1 

Hence, S k g = Pg^k £ £-9 to machine precision, and, to machine precision, 

Pglk+l = ^ E (t*, (7fc)j/6,i/ ) «7/6,I/a = ^7fc- (6-19) 

On the other hand the orthogonal component ^k^Pglk is per (6.16) multiplied 
by 2, i.e., to machine precision, 



7fc+i - P 9 lk + i = 2 (7* - ^7fc) • (6.20) 

As a consequence, the algorithm starts to diverge beyond the point where 
machine precision is reached. 

The observations just made continue to hold for the more general algorithms in 
Subsections 4.1 and 4.2. Thus no substantial after-convergence error build-up 
occurs for the algorithms of Subsection 4.1. For the algorithms of Subsection 
4.2, with basic recursion step 



m— 1 

7fc+i = E h mjZiik (6.21) 

3=0 

the terms with odd j all lie in C g , and those with even j are given within ma- 
chine precision by b m j^k- Since J2j even b m j = 2 m ~~ 1 , the out-of-space component 
in 7fc gets multiplied by 2 m ~ 1 in each iteration step. 
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7 Zak domain considerations 



We consider the case that ab =| with integer p, q > such that gcd (p, q) = 1, 
and we define the Zak transform Z as (the extension to L 2 (R) of) the mapping 



h^(Zh){t,u) = b~ 1 ' 2 £ h[ —^)e 2mku , t,veR. (7.1) 

fc=-oo \ " J 

We refer to [34] and to [17, Sec. 1.5], for more details on the Zak transform 
and its role in Gabor analysis. 

For f,h e L 2 (R) we set (when t, v e R) 



$ / (^)=p- 1/2 f(^)ft-^,^+-)) , (7.2) 

V V 9 P/ / fc=0,...,p-l,J=0,...,g-l 



and 



A f,h (t, v) = (t, '^)), ; . „ /; = $ / («, ^) (*, i/))* , (7.3) 

where the * denotes conjugate transpose. Now (g, a, b) is a Gabor frame, with 
frame bounds A > 0, < oo if and only if we have AI pxp < A" (t, v) < BI pxp 
for almost all t, v e R, with A and i? the largest and smallest positive real 
number for which the respective inequalities hold. The frame operator S of 
(g, a, b) is "represented" by A" through the formula 

$s/ = A93$ / ) f e L 2 (R) , (7.4) 

with matrix multiplication at each point (t, v) G R on the right-hand side 
of (7.4). This formula extends as follows. Assume that ip is continuous and 
positive on [A, B]. Then 

<&*W f = <p(A BB )&, feL 2 (R), (7.5) 

which is the basic formula for functional calculus in the Zak transform domain. 

We consider in this section the critical case a = b = 1 (in Sections 9 and 10 
more general rational ab will be dealt with). Then considerable simplifications 
occur since all the matrices $, A reduce to scalars. The formula (7.5) then 
becomes 
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(Z (tp (S) /)) (t, v)=<p (\(Zg) (t, u)\ 2 ) (Zf) (t, v), t,ueR, (7.6) 
for /eL 2 (R). In particular we have 



t,ueR, (7.7) 



(7.8) 

To illustrate the relevance for the algorithms, we consider algorithms II and 
IV for all scaling strategies. As to initial scaling, we assume that g and S are 
scaled such that (g, a — 1, b — 1) has best upper frame bound B < 2, which 
means that \Zg\ < 2 everywhere. We let 

G = Zg,Y k = Z lk . (7.9) 

Then by functional calculus in the Zak transform domain, the algorithms II 
and IV (initial scaling) assume the form 

r = G ; T k+1 = \^u-\ \Y k \ 2 T fc , k = 0, 1, ... , (7.10) 

and 



M^) = (z(^))M = fflM, 



(Zg d ) (t,u) = {z (S^gyj^u) 



\{Zg) (t,u) 
(Zg) (Mj 
\{Zg) 



(ZgT(t,uy 



t, v e 



r = g ; r fe+1 = 2r fe - |r fc | 2 c, fc = o, 1, . . . , (7.11) 

respectively, where the relations in (7.10) and (7.11) are to be considered at 
each point (t, v) G R. These recursions are then quite easily analyzed by ele- 
mentary means. For instance, one sees that the assumption B < 3 is necessary 
and sufficient for (7.10) to converge to exp(iarg(G)) everywhere, while the 
assumption B < 2 is necessary and sufficient for (7.11) to converge to 1/G* 
everywhere. Unbounded recursions result when we would have allowed B to 
be larger than 5 and 2, respectively. 

Next we consider algorithms II, IV using norm scaling so that (7.10) and (7.11) 
are to be replaced by 



r = G;r fc+1 = ^^-iH^Jl ? fc = 0,l,..., (7.12) 
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and 



r = G;r fc+1 = 2— , fc = 0,l,.... (7.13) 

II 1 fell IIUHI 

The norms used here are L 2 (JO, l) 2 )-norms. We consider the case that 



Zg = 1 on N ,Zg = x > on M, (7.14) 

where N, M are two measurable sets C [0, l) 2 such that N n M = 0, N U 
M = [0, l) 2 . Then (g,a— 1, b — 1) is a Gabor frame with best frame bounds 
A = min (1, x 2 ), B = max (1, x 2 ). Furthermore, 



Zg 1 = 1 on [0, l) 2 ; Zg d = - on M. (7.15) 

Ob 

We have for both algorithms II and IV that 



r fc = c fc on N , T k = d k on M, (7.16) 

where c k , 4 follow recursions that can be made completely explicit (using 
that, for instance, ||r fc || = ((1 — e) c 2 k + ed 2 k ) 1 ^ 2 , where e — /i (M)). Due to the 
norming operations in the recursion steps, either recursion stays bounded. 

We consider the case that e — /j, (M) <C 1. Then an elementary analysis shows 
the following: There is a 5 > such that for recursion (7.12), (7.16) there 
holds 



• x G (0, - 5 J => c k , 4 — > 1, 

ue + 5, - =>• c k — > 1, d k — > -1 

• x G (\/5 + 5, oo) =>■ chaotic behaviour. 

Also, there is a 5 > such that for recursion (7.13), (7.16) there holds 

• x e (0,V2- $) ^ c k ^ l,d k ^ ±, 

• x G (y/2 + 5, oo) =>• 4 < 0. 



8 Sampling and periodization of Gabor frames 



The algorithms considered in this paper and in [20,18,19] have been formulated 
for time-continuous Gabor frames while the tests we perform must take place in 
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a finite setting. The transition from continuous to discrete/finite Gabor frames 
by sampling and periodization has been discussed in [16] and later in [22,30], 
see also [7, Subsec. 8.4] and [6, Sees. 10.2 and 10.3]. Let a, M be positive 
integers, and assume that (g,a,l/M) is a Gabor frame with frame bounds 
A > 0, B < oo. Furthermore, assume that g satisfies the aforementioned 
condition A and the so-called condition R: 



R 

i r/ 2 

]=-oo 



oo rs/2 

lim ^ -/ \g{j + u) - g(j)\ 2 du = 0. 



The conditions R and A are not very restrictive; they are, for instance, satisfied 
by all members g of Feichtinger's algebra S , see [7, comment after Thm. 8.4.2]. 
Then the system 



9na, m /M ■= (gna,m/M ( r )) rgZ , U E Z,m = 0, . . . , M - 1 (8.2) 

is a discrete Gabor frame with frame bounds A, B, the dual window g d of the 
frame (g,a,l/M) satisfies conditions R and A, and the dual window (fi ,D ) 
corresponding to the discrete Gabor system in (8.2) is obtained by sampling 
g d : 



(9 D Y (r) = (g d ) D (r) , reZ. (8.3) 

The transition from discrete Gabor frames to discrete, periodic Gabor frames 
is just as convenient. Assume that we have a g e I 1 (Z) such that the discrete 



Gabor system \q n am/M) is a discrete Gabor frame with frame 

V y ' 1 JneZ,m=0,...,M-l 

bounds A > B < oo. Let L = Na = Mb for some positive integers N, b, and 
define 



9 P {r)= >T g(r-jL), r e Z. (8.4) 

3=-oo 

Then the system 

(C/MW) reZ - n = 0,...,iV-l,m = 0,...,M-l, (8.5) 

is a discrete, periodic Gabor system with frame bounds A, B, the dual win- 
dow g d of the discrete Gabor system is in I 1 (Z), and the dual window (g P ) 
corresponding to the discrete periodic Gabor system in (8.5) is obtained by 
periodizing g d : 
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d P 

(g P ) (r) = (g d ) (r) = J2 9 d (r- jL) , r G Z. 



(8.6) 



j=-oo 



An important extension of these results is given in [7, Subsec. 8.4]. Assume 
that ip is analytic in an open neighbourhood containing [A, B], where A > 0, 
B < oo are frame bounds of the Gabor frame (g,a,l/M) with g satisfying 
condition R and A. Then (p (S) g satisfies R and A as well, and 



where S D is the frame operator corresponding to the system in (8.2). The 
approach in [7, Subsec. 8.4] (which uses the Dunford representation of op- 
erators as well as theorems of the Wiener l//-type) can be mimicked so as 
to generalize the transition result from discrete Gabor systems as above with 
g G I 1 (Z) to discrete, periodic Gabor systems. Thus, with <p as above and 
\9nam/M) a discrete Gabor system with g G I 1 (Z), frame bounds 

V ' / n&,m=0,...,M— 1 

A > 0, B < oo and frame operator S, we have ip (S) g G I 1 (Z) and 



where S p is the frame operator of the system in (8.5). In particular, we see 
that the sampling-and-periodization approach is valid for the tight window g l 
in which case we should consider ip (s) = s -1 / 2 . 

It follows from the above results that the algorithms can be considered for 
discrete and for discrete, periodic Gabor frames. The findings for these systems 
are of direct relevance to the algorithms we have considered for the time- 
continuous case. 



9 Implementational aspects 

All implementations are done in the finite, discrete setting of Gabor frames. 
We denote for g G C L and a, b G N by (g, a, b) the collection of time-frequency 
shifted windows 




(8.7) 




(8.8) 



gna,mb, n<EZ,m<EZ, 



(9.1) 



where for j, k G Z we denote 



9j,k = e 



■ 2mkl/L g(i-j), J = o,...l-i. 



(9.2) 
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Note that it must hold that L = Na = Mb for some M, N G N. Additionally, 
we define c, d, p, q G N by 

c = gcd(a,M), d = gcd(b,N), V=~ = K Q= — = ^- (9-3) 

ca c d 

With these numbers, the density of the Gabor system can be written as(a&) / L = 
p/q, where p/q is a irreducible fraction. It holds that L = cdpq. 



9.1 Matrix representation and the SVD 



Let O g G C LxMN be the matrix representation of the synthesis operator of a 
Gabor frame so that 



( 9 )l, m+ nM = 9rna,nb(l), I = 0, . . . , L ~ 1, (9.4) 

for m — 0, . . . , M—l, n = 0, . . . , N—l. Hence O g has the column vectors g ma ,nb- 
The matrix representation of the frame operator corresponding to (g, a, b) is 
then given as O g O*. Since (g,a,b) is a frame we have that O g has full rank 
L < MN. 

Assume that ip is continuous an positive on a (S). From 



na,mb 

= (p(S) (9.5) 

we have that 

O ip(s)g = V (S)O g = V (O g O;)O g . (9.6) 

1 1 2 

Furthermore, note that for the Frobenius norm ||O ff || fro we have ||O ff || fro = 

1 1 1 1 2 

MN \\g\\ , since all columns of O g have norm \\g\\. 

The iterations can be written in terms of the synthesis operator matrices as 
follows. Denote the synthesis operator matrix corresponding to the Gabor 
frame (7^,0,6) by Q k . Then we can write the iteration step for algorithm II 
with norm scaling as 

^ = O g ;" fc+ i = -l?7^-- ,,i"n"?n"n » ^ = 0, 1, . . . . (9.7) 

We shall consider the thin SVD of the synthesis operator matrices. Thus we 
let O g = UYlV*, where U € C LxL is unitary (O p has full rank), E G M LxL 
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is a diagonal matrix with positive diagonal elements and V G C has 
orthonormal columns. With ip as above, we compute the thin SVD of 0^s) g 
as 

0^s )9 = V (O g O* g ) O g = p (UZ 2 U*) UEV* 

= Up(Y?) U*UYV* = U^p(Y?) V*. (9.8) 

Here we have used that p {UT?U*) — Up (£ 2 ) U*, a basic fact in the functional 
calculus of matrices. The equation (9.8) shows that 0^ S ) g has the same right 
and left singular vectors as O g , and the singular values transform according 
to ip — > up 2 (er). As a consequence we have, 

O gt = UV* , O gd = ITE^V*, (9.9) 

for the synthesis operators corresponding to ((?*, a, b) and (^g d , a, bj , respec- 
tively, for which we should take p> (s) = s -1 / 2 and p (s) = s^ 1 in (9.8). We 
thus see that we have obtained the matrices occurring in the polar decompo- 
sition of O g and the Moore-Penrose pseudo-inverse of O g . 

A further observation is that ||0 9 ||f ro = X^=i crj, where <7j, j — 0, . . . , L — 1, 
are the singular values of O g . Letting <pk,j be the singular values of Qk and 
using that 

n k = un k v\ 

we can write the iteration step in (9.7) on the level of singular values as 

3 1 1 1 

V ^j=i CT fc,j V ^i =1 

where j = 0, . . . , L — 1, and k — 0, 1, 

9.2 Factorization of finite, discrete Gabor systems 

Similar to the Zibulski-Zeevi representation of the Gabor frame operator in 
the continuous case, see [34,35,17], it is possible to compute the actions of 
the finite, discrete Gabor frame operator (and also the analysis and synthesis 
operators) very efficiently Several equivalent methods exists using almost the 
same number of operations, but differing in the order. In [4] a finite, discrete 
version of the Zibulski-Zeevi representation is developed. Another method was 
developed in [26] and [31]. Unfortunately, [31] contains some errors, which 
have been corrected in [3] . In the following we shall present the Zak-transform 
method from [4]. 
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For h G <C L and K G {0, ...,L- 1} such that £ G N, we define the finite, 
discrete Zak transform by 



(Z K h)(r,s) = J- ]T ^(r-ZK)e™/ L , r, s G Z. (9.11) 

The finite, discrete Zak transform is quasi-periodic in its first variable and 
periodic in the second, 



(Z K h) (r + kK,s + l^j = e 2 ™ ksK ' L (Z K h) (r, s) , (9.12) 

see [15] for more details. The values {Z^K) (r, s) of a finite, discrete Zak- 
transform on the fundamental domain r = 0, . . . , K — 1, s — 0, . . . , K/L — 1 
can be calculated efficiently by K FFT's of length K/L. To obtain values 
outside the fundamental domain, the quasi-periodicity relation (9.12) can be 
used. 

We define the cd matrices $£ s of size p x q and the p x p matrices A^ T by 



Ks = \M ((Zaf) (r + kM,s + W)) fc=0 ,..., p - M= o,... >? -i , (9-13) 
where r — 0, . . . , c — 1, s = 0, . . . , d — 1 and 

Mi={Mi)^ = P.") 

With these definitions is holds that the frame operator S of (g, a, b) is "repre- 
sented" by A" through the formula 

$sf = A 99Qf^ feC L , (9.15) 
see [4]. 

With this efficient representation of the frame operator of a finite, discrete 
Gabor system, we may express the iterations schemes in the finite, discrete 
Zak domain. We let 

G = ®>,r k = $v>,A k = = r k (r k y . (9.16) 

By functional calculus, algorithm II in the finite, discrete Zak transform takes 
the form 
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r — G ; T k+ i — -T k — -A k T k , k — 0, 1, 

The expressions for the other iterations types are similar. 



(9.17) 



9.3 Other methods 

We have considered two other methods of computing the canonical tight win- 
dow utilizing the factorization (9.13). To calculate the factorization of the 
canonical tight window g f ,^ 9 , we use an eigenvalue decomposition of the 
factorization of the frame operator of (g,a,b): For each r = 0, ...,c — 1, 
s — 0, . . . , d — 1 compute U riS , D r s such that A g r g s = U r ^ s D r ^U*^, where U TyS is 
unitary and D r s is diagonal and set 

('-J),,: 2 u;, s n, s - (9.18) 

We shall refer to this method as the EIG method. The other method uses 
(9.9) applied to the matrices of the factorization: For each r = 0, . . . , c — 1, 
s — 0, . . . , d— 1 compute U TyS) D r s , V rjS such that = U r ^D r>s V* s , where U r>s 
is unitary, D r s is diagonal and V r ^ s has orthonormal columns. Then it follows 
from functional calculus in the Zak transform domain (pretty much as in (9.8); 
also see (7.5)) that 

®C = U ryS V* s . (9.19) 
We shall refer to this method as the SVD method. 

For computing the canonical dual window we have considered simply inverting 
the matrices of the factorization of the frame operator: 

fo9 d = (A99 \~ l (h9 
r,s y r,s J r,s' 

We shall refer to this as the INV method. 

9.4 Implementational costs 

The computation of <3> 9 needs to be done before the iteration step. It can be 
computed using 5L log 2 d flops. This transforms the initial window g into the 
finite, discrete Zak domain. All computations in this domain are then done by 
multiplication of p x q and p x p matrices. The transform $ is unitary from 
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Table 3 



Method: 


Flop count per iteration: 


I. 


16Lp + |cdp 3 . 


II. 


16Lp. 


III. 


24Lp. 


IV. 


16Lp. 


V. 


24Lp + 8cdp 3 . 




Total flop count: 


INV. 


16Lp + §cdp 3 . 


EIG. 


24Lp + 14cdp 3 . 


SVD. 


64Lp + 32cdp 3 . 



This table shows the flop count of each of the considered methods. The flop count 
does not include the cost of the pre- and post-factorization. The application of 
an inverse matrix needed for the algorithms I and INV is done using a Cholesky 
factorization followed by two substitutions. An iteration step of V takes more 
flops to compute than an iteration step of III, because we need to compute the 
two terms S^g and Sk^k- The flop counts for EIG and SVD methods are only 
approximations, because eigenvalues and singular values can be calculated by 
many different methods with different flop counts, and because the process usually 
involves an iterative step, see [10]. 

C L with Euclidean norm into £. cxdxpxq , also with Euclidean norm. This gives 
an easy way to calculate the norms needed for the norm scaling. 

We count the number of real floating point operation needed, and assume 
that everything is done using complex arithmetics. The flop count for a single 
iteration step in the transform domain for each of the 5 algorithms can be seen 
in Table 3. 

A quick comparison show that the iterative methods for computing the tight 
window are comparable in number of flops to the EIG and SVD methods, if 
the number of necessary iterations is not to big. For the inverse iterations, 
the situation is different: Computing the inverse of the block matrices by a 
direct approach requires only slightly more flops than a single iteration step 
of algorithm IV, so an iterative method will always use more flops than the 
direct approach. However, there might be situations were it is not desirable to 
compute the inverse. For instance, if the initial window g has small support 
then the iteration steps can be performed by multiple passes through a filter 
bank. 



27 



9.5 Stopping criterion 



Because of the guaranteed quadratic/cubic convergence of the algorithms, it 
is possible to devise a simple yet powerful stopping criterion: We consider the 
difference 

ll7 " +1 ~7 fcl1 . (9.20) 
Il7fc+i|| 

When this difference is close to the machine precision eps, the algorithm con- 
sidered has converged. This is a standard stopping criteria, but using it this 
way means that we have done exactly one iteration step too much. Therefore, 
we stop when (9.20) is less that ^/eps and ^Jeps for the algorithms having 
quadratic and cubic convergence, respectively. 



9. 6 Window functions 



As the basic window functions we shall use the Gaussian (p e L 2 (R) and the 
hyperbolic secant ip e L 2 (R) given by 



< P(t) = (\) V4e ~^ teR, (9.21) 
rj} (t) = ^|sech (irt) , f el. (9.22) 

Both functions are invariant with respect to Fourier transformation. In order 
to generate a range of functions, we introduce a parameter w that dilates the 
functions by the unitary operator D w : L 2 (R) — > L 2 (R) given by 



(D w f)(t)=(wr 1 / 4 f(^=j, tel. (9.23) 
Applying this gives 



(D w <p) (t) = Vw (t) = (|) ~ 1/4 e"^ 2 / w , t e R, (9.24) 

T 1 ><-<h (t " 



(D w tP) (t) = i> w (t) = ^-w-^seeh ^t—J , tel. (9.25) 

It holds that the Fourier transform of ip w is (fii/ w and similarly for ip w . As 
window functions for the testing of the iterative algorithms we shall use finite, 
discrete versions of these, obtained by the sampling-and-periodization process 
described in Sec. 8: 
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Figure 2. The first row of graphs shows the Gaussian function (p± and the 
hyperbolic secant ip®. Both functions have length L = 120 and unit norm. The 
second row of graphs shows their canonical dual windows, and the third row shows 
their canonical tight windows. In the second column, the functions are shown on a 
logarithmic scale to display their decay. 




The properties from the continuous setting carry over: The functions have 
unit norm, and the Discrete Fourier Transform of </?^ is (fu w and similarly 
for ip®. For more details on the hyperbolic secant as a Gabor window, see 
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Figure 3. Fig. (a) shows the function MONSTER of length 600 with a single 
singular value set to cr rea i = 6. Fig. (b) shows the Discrete Fourier Transform of the 
function. 



[21]. Figure 2 shows a Gaussian and a hyperbolic secant and their respective 
canonical dual and tight windows. The different decay properties of the two 
bell-shaped functions are visible on a logarithmic scale. However, even though 
the Gaussian and the hyperbolic secant have different decay properties, their 
canonical windows have almost the same. 

To produce examples for which the norm scaling methods diverges, we have 
constructed a function (MONSTER) which is a Gaussian function modified in 
such a way that the first singular value, er rea i, of the matrix representation of the 
Gabor synthesis operator that corresponds to a real and symmetric singular 
vector, is given a large value. The function is shown on Figure 4. This function 
is a generalization to the case of rational oversampling of the counterexample 
given in Sec. 7 and exploits that the iterations can be considered as scalar 
iterations of the singular values of the Gabor synthesis operator, (9.10). 



10 Experiments 



This section contains the results from the experiments we have done in order to 
test the algorithms thoroughly and to demonstrate the various aspects of the 
algorithms that have been shown analytically. The computations have been 
done in Matlab and Octave, and the full source code is available for download 
from http : //www2 .mat . dtu. dk/people/P . Soendergaard/ iteralg/ We will 
show figures demonstrating the important aspects, but since we cannot include 
all material, the reader is encouraged to download the software and experi- 
ment. 
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(a) Tight iterations. 



(b) Dual iterations. 



Figure 4. The figure shows the behaviour of the 5 iteration types for the first 12 
iterations of each. The y-axis shows the Z 2 -norm of the difference between the 
iteration step and a precomputed, normalized solution. The system considered in 
the Gabor frame for C 432 , (</?f , 18, 18). The Gabor frame has a frame bound ratio 
of B/A = 2.03. 



10.1 Convergence and divergence of norm scaling 



Figure 4 show the convergence behaviour for a well-conditioned problem. The 
figure shows that algorithm I, II and IV exhibit quadratic convergence, III and 
V exhibit cubic convergence as proved in Subsections 4.1 and 4.2. Further- 
more, the algorithms for computing the tight window stay converged close to 
the machine precision, while the algorithms for computing the dual window 
diverge. Algorithm V also diverges faster than IV. This is as proved in Sub- 
section 6.3; the slopes of the two line segments beyond the 5 th iteration in 
Fig. 4(b) corresponds to divergence factors 2 (for IV) and 4 (for V). A visible 
numerical aspect is that iteration V is not able to reach full precision, because 
the iterand is quickly affected by the buildup of numerical errors. The conver- 
gence behaviour of the algorithms for the initial window being a hyperbolic 
secant is almost the same. 

Two examples of using different scaling strategies is show on Figure 5(a) and 
5(b). The figures show that initially scaling by the best scaling constant is as 
good as using norm scaling and using initial scaling by an easily computable 
scaling constant results in only 1-2 more iterations than using norm scaling. 
Comparing these methods to the method using optimal scaling, we see that for 
a well-conditioned problem then norm-scaling and optimal initial scaling are 
close to the optimal convergence. For a worse conditioned problem (Fig. 5(b)), 
optimal scaling clearly outperforms the other methods. However, this observa- 
tion has little practical relevance, because the computed canonical windows g d 
and g l will have a bad time-frequency localization. Higham [12] uses a scaling 
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Figure 5. Fig. (a) shows the convergence behaviour for algorithm II using four 
different scaling strategies: Norm scaling, initial scaling by the optimal constant, 
initial scaling by the dual lattice norm and constant optimal scaling. The system 
considered in the Gabor frame for C 432 , (</?f, 18, 18) . Fig (b) shows the same, but 
instead using the window <p^ 5 - This is a very narrow window, and the generated 
Gabor frame has a frame bound ratio of B/A = 180.8. 



strategy for algorithm I that approximates the optimal scaling. This requires 
an estimate for the smallest eigenvalue of the matrix, but this is easy to obtain 
since the matrix is inverted as part of the iteration step. For algorithms II- V 
we cannot use inversions, and so an estimate for the smallest eigenvalue (or 
lower frame bound) is difficult to obtain. We have therefore not pursued such 
a method for algorithms II- V. 



The iterations for computing the tight window are very robust when using 
norm scaling. It is easy to create examples of Gabor systems with frame bound 
ratios B/A > 10 12 for which the iterations converge, by using badly dilated 
Gaussians or by using a constant function with a small amount of noise added. 
However, by using the MONSTER function, it is possible to create an ex- 
ample for which the norm scaling iterations diverge. The behaviour of the 
dual lattice norm and ||<7 — 7fc|| 2 m each iteration step for a run of algorithm 
II step is shown in Fig. 6(a). It can be seen that the iteration converges to 
the wrong tight window. Another typical behaviour is that the iteration os- 
cillates between two different functions with the same dual lattice norm. The 
behaviour of algorithm IV on the same examples is shown in Fig. 6(b), the 
figure displays the optimal frame bounds of for each iteration step. Here 
we see exponential convergence of the lower frame bound of Zk to zero. 
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Figure 6. Fig. (a) shows the behaviour of the dual lattice norm and the P-norm of 
the difference between the iteration step and the normalized solution for a run of 
algorithm II using norm scaling. Fig. (b) shows the behaviour of the best upper 
and lower frame bound of in each iteration step for a run of algorithm IV using 
norm scaling. The system considered is in both cases are the Gabor frame for C 600 
with a = b = 20 using the MONSTER function. 
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Figure 7. The figure shows the numerical accuracy of three different methods to 
compute the canonical tight window. The plot is parametric in w. For each method, 
one line corresponds to a narrow window, w < 1, the other corresponds to a wide 
window, w > 1. Almost overlapping points on different lines corresponds to values 
w\,W2 such that w\ = l/w2- The error measure used is the dual lattice norm. The 
Gabor frames used are the Gabor frames for C 432 , ((f^, 18, 18). 
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10.2 Comparison with other methods 



Figure 7 show a comparison of the numerical precision of the algorithms for 
computing the canonical tight window compared to the numerical precision 
of other standard methods. The stability of the tight iterations proved in 
Subsection 6.3 is clearly visible. The method based on computing eigenvalues 
deteriorates quickly as the frame bound ratio increases, while the SVD be- 
haves much better. The eigenvalue method should not be used if the frame 
bound ratio of the problem is unknown. An explanation for this is that in the 
SVD method, the singular values are never considered, they are simple set to 
1. Therefore, roundoff errors on the small singular values do not affect the 
computation, in contrast to the EIG method, where round-off errors on the 
smallest eigenvalues are magnified because of the inversion of eigenvalues. For 
more details on the stability of computing eigenvalues and singular values see 

[!]■ 

The actual running time of the methods is determined by the flop count for 
each method (see the previous section for details) and of how fast the floating 
point operations can be executed by a computer. We will not give timings of 
the iterative algorithms, because we have not created optimal implementations 
of the algorithms, so timing them makes little sense. We note, however, that 
the key ingredients in the algorithms are FFTs of small length and matrix 
multiplications of small size matrices. Fast implementations exists for both al- 
gorithms, see [33,9]. This makes it possible to create efficient implementations 
of the iterative algorithms. 



10.3 Number of iterations 

To study how the number of necessary iterations depends on the frame bound 
ratio of the initial Gabor frame, we have plotted the number of iterations 
for the algorithms to converge, as a function of the frame bound ratio of the 
Gabor frame. Figure 8 shows such a plot, using Gaussians to generate Gabor 
frames with varying frame bound ratios. The jumps in the curves occur when, 
according to the stopping criteria, an additional iteration step is necessary. 
Even though algorithm III has cubic convergence, it is almost never able to 
compete with I. The jump in the graph for V is due to the magnification 
of round-off errors dominating the convergence, and causing divergence. The 
same happens for algorithm IV, but for considerable worse frame bound ratios 
(not visible on the graph). The graphs for the hyperbolic secant look similar. 

The number of necessary iterations might also depend on the size of the matrix 
blocks appearing in the factorization. This issue is slightly problematic to 
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Figure 8. The figure shows the number of iterations the algorithms need in order to 
reach machine precision. The three algorithms compared are I using norm scaling, 
and II-V using initial scaling by (6.9). The Gabor frames used are the Gabor 
frames for C 432 , (^,18,18). 



address, since creating a test problem involving bigger matrices also means 
altering the frame bound ratio. To minimize this effect, we have considered 
p, q running through the Fibonacci numbers, 2,3,5,8,..., such that p/q — > 
(y/E — l) /2. This creates a series of irreducible fractions p/q while keeping 
p/q close to a certain number away from 1. The result of the test is that 
the number of necessary iterations seems to be completely independent of the 
size of the matrix blocks! We have omitted the graphs, as they are simply 
horizontal lines. For algorithm I, it is proved in [23] that this is indeed the 
case. 



10. 4 Choosing an initial scaling 



Figure 1 shows that for each iteration type there is a range of values of the 
upper frame bound of the scaled window, B sca i ed , that will guarantee conver- 
gence. Figure 9 shows an example of the effect of prescaling the input window 
to obtain specific values of B sca i ed . As shown on Fig. 1, algorithm III diverges 
if B SC aied is larger than 7/3. The dual iterations IV and V diverge if B sca i ed > 2 
and II has a chaotic behaviour if 3 < B sca i ed < 5 and diverges if B sca i ed > 5 
(not shown on the plot). 

The choice of B that minimizes the number of iterations is to choose B ac- 
cording to Table 1. An estimate for this is difficult to calculate, as it involves 
an estimate for the lower frame bound. Fortunately, as can be seen on Fig. 
9 there is a large region around the optimal scaling point, where only 1 or 2 
extra iterations are needed. 
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Figure 9. The figure shows the number of necessary iterations to find the tight /dual 
window of a Gabor frame for C 432 , (y?f\ 18, 18) , as a function of the best upper 
frame bounds of the initial window. 

A A result on Condition A' 



Proposition 2 Assume that (g,a,b= 1/M) is a Gabor frame that satisfies 
condition A, where a, M e N. Also assume that <p is analytic around [A, B] 
and positive on [A, B], where A > 0, B < oo are lower, upper frame bounds of 
(g, a, b). Finally, let 7 = ip (S) g where S is the frame operator corresponding 
to (g,a,b). Then #,7 satisfies the Condition A', i.e., 



Yl I {9 , lj/b,l/a) 



hi 



< OO. 



(A.l) 



PROOF. We have for f,heL 2 



that 



(/, lna,mb) (g na ,mb, h) = ^ (/, (if (S) g) namb ) (g n a,mb, h) 
n,m n,m 

= J2(V ( S ) /> 9na, mb ) (gna,mb, h) = (Sip (S) f(h)2) 



We know that (7, a, 6)is a Gabor frame. Now take f,h G L 2 (R) such that 
(/, a, b) and (h, a, b) have finite upper frame bounds. Then by the fundamental 
identity of Gabor analysis, see [17, Subsecs. 1.4.1 and 1.4.2], 



(/) lna,mb) {9na,mb> tl) — (j) , ^j/b,l/a ) (fj/b,l/a, h) , 

3 J 



(A.3) 
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with absolute convergence on either side of (A. 3) . Thus Sip (S) has the dual 
lattice representation 



Sf(S) = ^(gn j m/a)u jJl . (A.4) 

a ° 3,1 

Now let ip (s) = sip (s). This ip is analytic around [A, B] and positive on [A, B]. 
By functional calculus of frame operators in the time-frequency domain, see 
[7, Sec. 8.3], there holds that Sip (S) has also the dual lattice representation 



Sip(S)^U(-HH*)) U hh (A.5) 

where H is the analysis operator with respect to the dual lattice, defined for 
f E L 2 (JR.), by 



Hf=((f, 9j/b , l/a ))^. (A.6) 

It follows from the proof of [7, Thm. 4.3], in particular from uniform bound- 
edness of (8.4.14) (with ip instead of ip), that 



E 



ab 



-HH* 



0,0;j, 



< oo. 



(A.7) 



By uniqueness of the coefficients in the dual lattice representation (just con- 
sider a well-behaved h such that {h na , m } ) ) n mgZ is a tight frame, i.e. such that 

Ujjh, j, I G Z, is an orthogonal set of functions), it follows that (fiS lj/b,i/a) 
oo, as required. 



< 



Note 1 It is implicit in the statement and proof of [7, Thm 4-3] that the 7 of 
the above result is such that condition A is satisfied by (7,0,6). 
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